suppressPackageStartupMessages({
    library(Seurat)
    library(cowplot)
    library(ggplot2)
})

library(future)
plan("multiprocess", workers = 4)
options(future.globals.maxSize = 24000 * 1024^2)

# remove.packages("Matrix")
# packageurl <- "https://cran.r-project.org/src/contrib/Archive/Matrix/Matrix_1.3-2.tar.gz"
# install.packages(packageurl, repos=NULL, type="source")


alldata <- readRDS("../analysis/filtered_IMM_DN.rds")

DefaultAssay(alldata) = "RNA"

1 Find variable genes

suppressWarnings(suppressMessages(alldata <- FindVariableFeatures(alldata, selection.method = "vst", 
    nfeatures = 2000, verbose = FALSE, assay = "RNA")))
top20 <- head(VariableFeatures(alldata), 20)

LabelPoints(plot = VariableFeaturePlot(alldata), points = top20, repel = TRUE)

2 Z-score transformation

alldata <- ScaleData(alldata, vars.to.regress = c("percent_mito", "nFeature_RNA"), 
    assay = "RNA")

3 PCA

alldata <- RunPCA(alldata, npcs = 50, verbose = F)

plot_grid(ncol = 3, 
          DimPlot(alldata, reduction = "pca", group.by = "orig.ident", dims = 1:2) + NoLegend(),
          DimPlot(alldata, reduction = "pca", group.by = "orig.ident", dims = 3:4) + NoLegend(),
          DimPlot(alldata, reduction = "pca", group.by = "orig.ident", dims = 5:6) + NoLegend())

VizDimLoadings(alldata, dims = 1:5, reduction = "pca", ncol = 5, balanced = T)

ElbowPlot(alldata, reduction = "pca", ndims = 50)

4 UMAP

5 Save data

saveRDS(alldata, file = "../analysis/filtered_IMM_DN_dr.rds")